N94-21475 

OPTIMAL CONVOLUTION SOR ACCELERATION 
OF WAVEFORM RELAXATION WITH APPLICATION 
TO SEMICONDUCTOR DEVICE SIMULATION 


5 // - a/ 


Mark Reichelt 

Research Laboratory of Electronics, Massachusetts Institute of Technology P / j 

Cambridge, MA < 


Summary 


In this paper we describe a novel generalized SOR algorithm for accelerating the 
convergence of the dynamic iteration method known as waveform relaxation. A new 
convolution SOR algorithm is presented, along with a theorem for determining the 
optimal convolution SOR parameter. Both analytic and experimental results are given 
to demonstrate that the convergence of the convolution SOR algorithm is substantially 
faster than that of the more obvious frequency-independent waveform SOR algorithm. 
Finally, to demonstrate the general applicability of this new method, it is used to solve 
the differential-algebraic system generated by spatial discretization of the time-dependent 
semiconductor device equations. 


Introduction 


To achieve highest performance on a parallel computer, a numerical algorithm must 
avoid frequent parallel synchronization [1], The waveform relaxation approach to solving 
time-dependent initial-value problems is just such a method, as the iterates are waveforms 
over an interval, rather than single timepoints [2, 3, 4]. Like any relaxation scheme, 
efficiency depends on rapid convergence, and there have been several investigations into how 
to accelerate WR [2, 5], including using multigrid [6] and conjugate direction techniques [7]. 

In this paper, we investigate using successive overrelaxation (SOR) to accelerate WR 
convergence. In particular, we show that the pessimistic results about waveform SOR 
derived in [2] can be substantially improved by replacing multiplication with a fixed SOR 
parameter by convolution with an SOR kernel. We derive the optimal SOR kernel using 
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Fourier analysis techniques and then demonstrate the effectiveness of the approach for a 
model parabolic problem. Finally, we demonstrate the general applicability of the approach 
by using the method to solve the time-dependent drift-diffusion equations associated with 
modeling semiconductor devices. 

We begin in Section 2 by reviewing waveform SOR, and in Section 3 we relate the 
algorithm to pointwise SOR to demonstrate the difficulty in accelerating WR with a fixed 
SOR parameter. In Section 4, we use Fourier analysis to derive the SOR kernel for the 
continuous WR algorithm, and give a proof of optimality. In Section 5 we briefly consider 
the effect of time-discretization, and in Section 6 we apply the method to device simulation. 
Finally, conclusions and acknowledgements axe given in Section 7. 


Waveform SOR 


In this section, we consider applying waveform relaxation methods to the model linear 
initial-value problem 

( 1 ) (ft + A ) *(0 — b ( t ) *(0) = ae 0 . 

where A E ! nxn , b(t) G E n is a given time-dependent right-hand side vector, x(t) E ffi n is 
the unknown vector to be computed over simulation interval t E [0, T], and xq € E 7 * is an 
initial condition. 

Given the relaxation splitting A = D — L — U, and subtracting successive waveform 
relaxation iterations, the waveform Gauss-Jacobi (WGJ) and waveform Gauss-Seidel 
(WGS) iteration equations, respectively, may be written as: 

(2) (£+D)Ax m (t) = (L + U) Ax k (t) 

(3) (i+D- L) A **•*(«) = U A 

where A« fcfl {t) = x m (t) — x k (t) is used to eliminate the right hand side b(t). 

The waveform SOR method for acceleration of WGS is a simple extension of algebraic 
SOR. To derive the waveform SOR iteration equation, compute a waveform xf 44 (t) on 
*€[0,21, as in WGS: 

( 4 ) (ft +a»)£f fl (£) = bi(t) - it, a i 3 x ’t{ t ) with xf fl (0)=x Q , 

3 = 1 j=i+i 

and then update x k {t) in the iteration direction by multiplication with an overrelaxation 
parameter u>, 

( 5 ) (*)«“*<(*)+«• [Zi* (*) - x i (*)] • 




• f -. • i r. 
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Combining equations (4) and (5) yields 


( 6 ) 


(f + °«) _ 

(1 - w ) [ (ft + a ii) xk i (*)] + U 


i - 1 


bi{t) - Y, a n x1 t (*) - J2 a n x i (*) 

3 = 1 j=»+l 


which, after subtracting successive waveform relaxation iterations, leads to 

(7) (i +D- uL) A x m (t) = [(1 - u){f t + D)+u>U] A x k (t), 
where Ajc** 1 ^) = x M (t) — x k (t). 

Note that the iteration matrices implied by equations (2), (3) and (7) correspond 
exactly to the standard algebraic relaxation and SOR matrices with diagonal matrix D 
replaced by 4- D ). Also note that waveform SOR as defined by (7) is not the same as 
the dynamic SOR iteration considered in [2], because, unlike WGJ or WGS, the waveform 
SOR iteration equations are not of the form 

(8) f t Ax k+1 + MAx k+1 = NAx k 


where M,N € M nxn . 


Relation to Pointwise SOR 


Discretizing (1) in time using a multistep integration method yields 

(9) ajxlm - j ] = ( b[m - j] - Ax[m - j]) , 

j = o 

where ao = 1 and x[m] denotes x(t) at timepoint t = mh with timestep h. Thus, the 
time-discretized model problem can be rewritten as a sequence of linear algebraic problems 

[/ + h/3 0 A] x [m] = 

s s 

(10) hp 0 b\m] - ^2 otjx\m - j] + h Pj (M m “ j\ ~ Ax[m - j ]) . 

i=i j=i 

We now compare the convergence of the waveform SOR method to the convergence of 
pointwise SOR, in which algebraic SOR is used to solve the matrix problem at each 
timepoint. 

The pointwise SOR iteration equations are derived by applying the relaxation splitting 
A = D - L - U to equation (10) and taking the difference between the (A:+l)st and fcth 
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iterations. More precisely, the pointwise SOR iteration equation applied to solve (10) for 
= a^fm] — x k [m ] is 

[ (I + hfoD) - uhfoL] A* w [ro] = 

(11) [ (1 - u) (I + hfoD ) + uhfoU] 

where u is the SOR parameter. It follows that the spectral radius of the iteration matrix 
generated by pointwise SOR at the mth timestep is 

(12) p ([ (I + hfoD) - uhfoL] _1 [ (1 - u j) (I + hfoD) + uhfoufj . 

If waveform SOR is used to solve the model problem (1), and a multistep method is 
used to solve iteration equation (7), then A x M [m], now denoting the discretized difference 
in waveform iterates, satisfies 

[Aa^fm - j } - (1 - u) Ax k [m - ;]] = 

j=0 

(13) ujL) A x M [m - j] + [(1 - w) D + u>U) A* fc [m -j]}. 
i= o 

This can be rewritten as the discrete-time analogue of (7): 

8 

y [ ( ctjl 4- hfijD) — uhPjL\Ax h¥1 \m — j] = 

3 = 0 

8 

(14) y [ (1 - u) {ctjl + hfljD) + uhPjU] Ax k [m - j]. 

3=0 

As the similarities of equations (11) and (14) suggest, if the time interval is finite, 
i.e. the number of timesteps is some finite L, then for a given timestep h and a given 
SOR parameter cu, the time-discretized waveform SOR method has the same asymptotic 
convergence rate as pointwise SOR. 

Theorem 3.1. On a finite simulation interval, the iterations defined by (11) and (14) 
have the same asymptotic convergence rate. 

Proof. Let y k denote the large vector consisting of the concatenation of vectors A x k [m] 
at all L discrete timepoints, i.e. y k = |Aa? fc [l] T , . . . , A« fc [L] T J . Collecting together the 
equations (14) generated at each timepoint into one large matrix equation in terms of 
vectors y k+1 and y k yields MAy k+1 = NAy k where M,7V 6 l LnxLn are block lower 
triangular banded matrices, with blocks of size n x n, and with block bandwidth s. It is 
then easily seen that M~ 1 N is block lower triangular, with diagonal blocks equal to 

(15) [ (/ + hfoD) - uh/3 0 L] ^ [ (1 - u ) (I + hfoD) + uhfou] . 
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Therefore, p(M l N) is given by (12), implying that the iterations defined by (11) and (14) 
have identical asymptotic convergence rates. C 


Theorem 3.1 suggests that parameter u> for waveform SOR should be chosen to be 
precisely equal to the optimum parameter for the pointwise SOR method. However, this 
does not necessarily lead to fast convergence, as the following example illustrates. 

Example 3.1. Let t G [0,2048], *(0) = 0, and let matrix A G K 32 * 32 and time- 
dependent input vector b(t) G t 32 of the model problem (1) be given by 


(16) 


A = 


b(t) = 


2 -1 

-1 2 


0 

0 


-1 

-1 2 

, , M ( l-cos(^) if t < 256 

where bi(t) = < \256/ 

I 0 otherwise. 


Consider the four problems generated by discretizing in time with the first-order backward 
difference formula, using 64, 128, 256, and 512 uniform timesteps of size h = 32, 16, 8 and 
4 respectively. 


Since the tridiagonal matrix A is symmetric and is consistently ordered [8, 9], the 
matrix (I+hfioA) of the pointwise time-discretized model problem (10) is also consistently 
ordered, and the optimum pointwise SOR parameter ujgpt is given by 


(17) 


W opt — 



where p\ = p(H G j) is the spectral radius of the pointwise Gauss-Jacobi iteration matrix 
Hqj = (I + hPoD)~ l (hPoL + hfi^U). For the four problems with 64, 128, 256 and 
512 timesteps, the optimum pointwise parameters uj opt are 1.669, 1.586, 1.482 and 1.364 
respectively. 

Curves PT64, PT128, PT256 and PT512 of Figure 1 show the convergence of the 
waveform SOR method versus iteration for the four problems with their optimum pointwise 
SOR parameters Uop t . Note that as the total number of timesteps is increased, the initial 
convergence rate is slower, approaching a limiting value of the convergence rate of the 
continuous Gauss-Seidel WR algorithm (shown as WR in Figure 1). In each case, the 
convergence rate of the waveform SOR eventually approaches the expected asymptotic 
value of uj opt — 1. Note that with a reasonable error accuracy tolerance such as 10~ 6 as a 
stopping point, the asymptotic convergence rate is ne ver reached. For comparison, Figure 1 
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also shows the superposition of four convergence plots (CSOR) of the new convolution SOR 
method to be introduced in the following sections. 



Fig. 1. Convergence of waveform SOR using 
the pointwise optimal parameter (PT) compared to 
waveform relaxation (WR), and convolution SOR 
(CSOR), with 64 , 128, 256 and 512 timesteps. 



Fig. 2. Effect on convergence of the 256-timestep 
waveform SOR of varying the SOR parameter from 
the pointwise optimum L 0 opt — 1.482. 


To illustrate the effect of choosing a different SOR parameter u, Figure 2 shows the 
convergence versus iteration of the 256-timestep example for waveform SOR with values 
of the SOR parameter u> not equal to the pointwise optimum ujop t = 1.482. When 
u = 1.30 < <jJ 0 pt, the convergence curve lies between the pointwise optimum curve and 
the WR convergence curve, i.e. both initial and asymptotic convergence rates are slower. 
By increasing the SOR parameter to w = 1.63 > uj opt , the initial convergence rate can 
be made faster at the expense of slowing down the asymptotic convergence rate. But as 
the u) = 1.70 curve shows, once the SOR parameter is increased beyond some point, the 
waveform SOR method may appear to diverge before eventually converging. Also, the 
solution produced by the ui = 1.70 example contains spurious oscillations, as shown in 
Figure 3. Note both the growth and translation of the oscillation with iteration. 

The optimum pointwise SOR parameter uJ opt does not dramatically improve the 
convergence rate of waveform SOR because the matrix M~ l N which describes the 
waveform SOR convergence is far from normal. This suggests that although the spectral 
radius of the iteration matrix determines the asymptotic convergence rate of waveform 
SOR, it does not determine the practically observable convergence rate. The convergence 
rate could be characterized, for example, by computing the pseudo-eigenvalues [10] of the 
waveform SOR iteration matrix. In the following section, we take an alternate approach. 
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Fig. 3. Delta waveform = xfg 1 (t) - 

xf 6 (t) versus time after iterations 250 and 500 , 
for the 256-timestep waveform SOR method using 
cu = 1.70, showing the growth and translation of an 
oscillating region . 



Fig. 4. The spectral radii as functions of frequency 
fi of the Gauss-Jacobi WR (solid), Gauss-Seidel 
WR (dashed) and waveform SOR (dotted) iteration 
matrices for an 8 x 8 version of the continuous-time 
problem of Example 3.1. 


Fourier Analysis 


In [2], the spectral radius of dynamic iteration operators which map x k to x M , such 
as those given by equations (2), (3), and (8), was related to their Fourier transform. In this 
section, we make a more detailed use of Fourier analysis to derive a frequency-dependent 
SOR parameter for the waveform SOR operator of equation (7). 

The Fourier transform of x k (t) is given by 

(18) x k (iQ) = r x k {t)e- in dt = F{x k {t)}, 

J— OO 

where 0 is frequency. Standard Fourier identities can be used to show that = 

H(iQ) Ax k (iQ), where for WGJ (2), WGS (3) and waveform SOR (7), the iteration 
operator H (iQ) is given by 

(19) H G j{in ) = (iftl + D)-\L + U) 

(20) H GS (ity = (ini + D-L)~ l U 

(21) HsoRiiV) = (MI + D-uL)- 1 {(l-uj)(ini + D)+u>U} 


respectively. The obvious interpretation of equations (19)-(21) is that the spectral radius 
p(H(in)) yields the asymptotic convergence rate for errors in the frequency component Q. 
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Figure 4 is a plot of the spectral radii of H GJ (iQ.), H GS {ity and H S or{ *0) for an 
8x8 version of the continuous-time problem given in Example 3.1, using u> = 1.49 for 
Hsodin). From the plot it is clear that very high frequency components of the error are 
damped much more quickly than low frequency components. However, the spectral radius 
p(H sor(^)) is greater than one over a range of frequencies, and therefore the waveform 
SOR iteration magnifies errors in this frequency range. This effect was predicted in [2] and 
is easily seen in Figure 3. 

This situation can be remedied by using a generalized SOR algorithm, in which 
equation (5) is replaced by an overrelaxation convolution with a time-dependent SOR 
parameter u>(t), 

/ OO r -1 

u>(t - t) • [^(t)— x£(t)] dr. 

■OO L 

The Fourier transform of the SOR operator is then given by 

(23) H G (itt) = [ifiJ + D — u(iCl)Lj J(l — + £>) + 

where u>(zQ) is the Fourier transform of the time-dependent w(t). We refer to the SOR 
algorithm represented by iteration matrix (23) as the convolution SOR algorithm (CSOR). 
The theorem below, which is the main result of this paper, gives a formula for determining 
the optimal frequency-dependent SOR parameter 


Theorem 4-2. If the spectrum of H GJ (iQ, ) lies on the line segment [— //i(ifi),/ii(ift)] 
with |/zi| < 1, then the spectral radius of H c (iO) is minimized at frequency by a unique 
optimum u>(iQ) = u op t iiOl) € C given by 


(24) 


u op t(iO) 


2 

1 -1- ^1 -/ii(ifl) 2 


where yf- denotes the root with the positive real part. 

Proof. For brevity, the argument (iQ) will be omitted in the following, and H G (u>) will 
denote the convolution SOR operator (at frequency Q) computed using SOR parameter u>. 

Let /ij = rifii denote each eigenvalue of H GJ , where G [—1,1]. Classical SOR 
theory [8, 9] guarantees that for each /Zj = 7\pi, there is an eigenvalue A* of H g (uj) which 
satisfies 

(25) A i - urifr + (w - 1) = 0, 
and therefore, from the quadratic formula, 

(26) + 
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Let u to be the conjectured optimal u> op t ■ Combining equation (24) with (26) yields 

(27) ^/|Aj| = \p\^opt r i + ‘\J r i ~ 1 = » 

where the rightmost equality follows from the fact that Jr t + ijrf — l| = 1 for E [-1, 1]. 
And as (27) holds for all i, 

(28) p{H ci^opt)) = I Ai| = I opt) I = l^opt ~ 1| • 


Equation (28) implies that p(Hc(u)) cannot be decreased below p(H ci^opt)) by using 
an u such that |u; — 1| > \^opt ~ 1|- This follows from the fact that, in general [8, 9], 

(29) p(H c (u>)) > |w- 1| 


for any u. 

To show that p(H c {u>)) also cannot be decreased by choosing a value of u such that 
|u/ — 1| < | ujopt — 1|, consider the eigenvalue A j corresponding to p\\ 


(30) 


y/Xj = f+{u) = 



and note that /+ : C -* C, given by equation (30), is a single-valued, continuous function 
that is analytic except at 


(31) 


Ui,Ul 2 


2 


1 ± 



Since |^i| < 1, points u\ and u> 2 lie in the interior and exterior, respectively, of the 
circle \lj — 1| = 1 in the complex tu-plane. Note that u)\ equals the conjectured tOopt from 
equation (24). 

Let D denote the interior of the curve given by the perimeter of the circle \u — 1| = 1, 
except with a cut along the line defined by the circle’s center and ui\ . The cut follows the 
line from the perimeter down to uq, and then back up the other side to the perimeter, as 
shown in Figure 5. The function /+ is nonzero everywhere within D, since equation (25) 
implies that a zero can occur only at w — 1, and /+(1) = p-\- Therefore, the minimum 
modulus theorem [11] implies that |/ + (u>)| attains its minimum value somewhere on the 
boundary of D. Finally, the lower bound in (29) implies that uq = u) op t in (24) is the only 
point on D which can achieve as low a p{H c{u)) as given in (28), completing the proof. □ 


Note that when the eigenvalues p, he on a real fine segment, this is yet another 
alternative proof of a classic SOR Theorem [8, 9, 12]. Also note that, in general, the 
optimal overrelaxation parameter u is complex. 
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Fig. 5. The region D and branch cuts in the complex w-plane. 

The conditions of optimal SOR parameter Theorem 4.2 are satisfied by a large class of 
matrices. 


Corollary 4-3. If A in (1) is consistently ordered, symmetric, and has constant diagonal 
D = dl, then the optimal SOR parameter is given by (24), 

2 


(32) 


w 0 pt(ifi) = 


1 + 


\ 


_ ( dfj, i V 

\d> “h J 


where denotes the spectral radius of£> 1 (L + U). 


Proof. To show that Theorem 4.2 applies, note that (19) implies that for a constant 
diagonal A, the eigenvalues are given by = d/M)/(d + if)), where /x 0 

are the eigenvalues of D~ X (L + U ). Since fi o lie on the real axis, the /x(if)) lie on a line 
rotated in the complex plane. □ 


Corollary 4-4- If A in (1) is consistently ordered, symmetric, and has constant diagonal 
D = dl, then the optimal time-dependent SOR convolution waveform u{t) is real. 

Proof. Equation (32) implies that o>op t (if)) is a conjugate-symmetric function of f). □ 


Discrete-time Modification 


For the sake of brevity, we consider only the first-order backward difference formula, 
in which case equation (14) becomes 

4- h(D — wZl)Aa; fc+4 [m] — Ajc^fra — 1] = 

(33) (1 — u) Aa; fc [m] + [ (1 — u>) hD + hujU^A.x k [m] — (1 — u>) A x k [m — 1], 
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where h is the uniform timestep. The z-transform of x[m\, defined by 


(34) x{z)= Y x[rn}z m = Z{x[m}}, 

TV=~OQ 

may be used to show that A x m (z) = H c {z ) A x k {z), where the ^-dependent convolution 
SOR operator is 

H c {z) = - I + D- uj(z)L^j (l - w(z)) ^ + “( Z ) U • 

Since u>(z) depends on z, overrelaxation becomes a convolution sum 

o° 

(35) x? 1 [m)*-x k i {m}+ Y “i m ~ j ] ' W) » 

i=-oo 

where uj(z) — Z{u[m]}. To determine the optimal a ;(z), we have the following theorem, 
whose proof is analogous to that of Theorem 4.2. 

Theorem 5.5. If the spectrum of H GJ {z ) lies on the line segment [-/xi(z),M 2 )] 

|^i| < 1, then the spectral radius of H G {z ) is minimized at z by the unique optimum 
u(z) = u) op t(z) (EC given by 

(36) opt(z ) 
where yf denotes the root with the positive real part. 

In Example 3.1, matrix A has constant diagonal D = dl, so that 

(37) w opt(z ) 


where pi denotes the spectral radius of D~ l {L + U). Thus, to compute the optimal 
convolution SOR sequence u[m] for the four CSOR plots of Figure 1, equation (37) was 
used to compute u ;(z), and then the inverse z-transform of u(z) was computed analytically 
by series expansion. 
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Device Transient Simulation 


A device is assumed to be governed by the Poisson equation, and the electron and hole 
continuity equations: 

(38) V 2 u + c\{p — n + N) = 0 

dTh 

(39) V 2 n - VnVu - nV 2 u = c 2 — 

at 

(40) V 2 p + VpVu + pV 2 u = C3— 

at 

where u is the normalized electrostatic potential, n and p are the electron and hole 
concentrations, N is a background concentration, and 01,02,03 are physical constants [13]. 

Given a rectangular mesh that covers a two-dimensional slice of a MOSFET, a common 
approach to spatially discretizing the device equations is to use a finite-difference formula 
to discretize the Poisson equation, and an exponentially-fit finite-difference formula to 
discretize the continuity equations [13]. On an iV-node rectangular mesh, the spatial 
discretization yields a differential-algebraic system of 3 N equations in 3 N unknowns. 

The convolution SOR method was implemented in the WR-based device transient 
simulation program WORDS [14]. WORDS uses red/black block Gauss-Seidel WR, where 
the blocks correspond to vertical mesh fines. The equations governing nodes in the same 
block are solved simultaneously using the first order backward-difference formula. The 
implicit algebraic systems generated by the backward difference formula are solved with 
Newton’s method, and the linear equation systems generated by Newton’s method are 
solved with sparse Gaussian elimination. 

The three MOS devices of Figure 6 were used to construct six simulation examples, 
each device being subjected to either a drain voltage pulse with the gate held high (the D 
examples), or a gate voltage pulse with the drain held high (the G examples). All examples 
ranged from low to high drain current, and in the G examples, the gate displacement current 
was substantial because the applied voltage pulses changed at a rate of .2 ~ 2 volts per 
picosecond. 


device 

description 

mesh 

5 v 

unknowns *~\ 5 | V j — \ 

a 

kar 

ldd 

soi 

abrupt junction 
lightly-doped drain 
silicon-on-insulator 

19 x 31 
15 x 20 
18 x 24 

1379 

656 

856 

1 1 

< > 

2.2 microns 

“T 

0 psec 


Fig. 6. Description of devices and illustration of the drain-driven karD example. 
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Figure 7 shows the convergence of the six examples as a function of iteration for WR, 
ordinary waveform SOR (using the pointwise optimum parameter), and the convolution 
SOR algorithm. The convolution SOR sequence u)[m] was calculated by linearizing (38)- 
(40) about the initial condition, estimating the spectral radius of the iteration matrix as a 
function of z, applying Theorem 5.5 and inverse transforming. Both overrelaxation methods 
were applied only to the potential variable u. All simulations began with 64 initial WR 
iterations, and used 256 equally-spaced timesteps. In Figure 7, convergence was measured 
using the terminal current error. 

Despite the nonlinearity of the semiconductor equations, the convolution SOR algo- 
rithm converged substantially faster than either WR or ordinary waveform SOR, demon- 
strating the robustness of the approach. 




0 100 200 300 

iterations 


0 100 200 300 

iterations 





Fig. 7. Terminal current error of the six examples as a function of iteration for WR (dashed), ordinary 
waveform SOR (dotted), and convolution SOR (solid). 








Conclusion 


In this paper, a new waveform overrelaxation algorithm was presented and applied 
to solving the differential-algebraic system generated by spatial discretization of the time- 
dependent semiconductor device equations. In the experiments included, the convolution 
SOR algorithm converged robustly, and substantially faster than ordinaiy WR. 

The author would like to acknowledge extensive conversations with his advisor, 
Professor Jacob White, and also thank Professors Alar Toomre, Donald Rose, Paul 
Lanzcron, Andrew Lumsdaine and Olavi Nevanlinna for many valuable suggestions. 
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